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An  Accurate  Method  for  FD-TD  Plane  Wave  Propagation 
for  Wide  Band  Excitations 


J.  T.  MacGillivray  and  D.  Dietz 
Air  Force  Research  Laboratory,  Directed  Energy  Directorate 
Kirtland  AFB,  NM  87117-5776 

Abstract 

In  this  paper  an  accurate  and  efficient  technique  is  presented  for  propagating  wide  band  plane  waves  in  a  three- 
dimensional  (3-D)  finite  difference  time  domain  (FD-TD)  grid.  Commonly,  the  incident  pulse  is  created  via  a 
separate,  one  dimensional  (1-D)  auxiliary  field  grid  and  then  applied  to  the  3-D  grid.  Although  this  paper  focuses 
on  the  separate  field  formalism,  the  method  can  be  applied  to  the  scattered  field  formalism  as  well.  Ai  the  wave  is 
applied,  a  3-D  numerical  wave  is  excited  and  propagates.  If  the  wave  number  band  of  the  auxiliary  field  grid  is  not 
nearly  equal  to  that  of  the  3-D  numerical  wave,  error  accumulates  resulting  in  inaccurate  results.  However,  by 
reconstructing  the  auxiliary  field  grid,  its  wave  number  band  can  be  closely  matched  to  that  of  the  3-D  numerical 
wave  yielding  a  much  higher  fidelity  propagating  plane  wave. 


I,  INTRODUCTION 

With  the  advent  of  much  faster  computer  CPU’s,  larger  processor  main  memory  and  massively  parallel 
machines,  it  is  now  possible  to  computationally  model  electrically  large  media  such  as  building  environments,  with 
techniques  that  previously  were  not  practical.  Understanding  the  microwave  environment  inside  a  building  is  of 
much  interest  for  communications,  for  example.  Features  such  as  reinforcement  bars,  I-beams,  conduit  and  interior 
equipment  scatter  fields  creating  a  challenging  task  for  an  analyst  to  assimilate  and  produce  useful  information.  It  is 
in  this  arena  that  time  domain  computer  modeling  techniques  can  contribute  significantly  by  calculating  overall  field 
distributions  and  areas  of  high  and  low  energy.  Arrays  of  virtual  field  sensors  record  data  that  can  be  used  to 
generate  animated  graphics  of  electromagnetic  fields;  this  is  not  possible  in  general  with  laboratory  measurement 
data.  Recent  research  into  building  microwave  environments  at  the  Air  Force  Research  Lab  (AFRL)  has  driven  the 
need  to  model  plane  wave  propagation  of  wide  band  pulses  through  electrically  large  meshes.  This  need  has 
stimulated  modeling  research  to  investigate  numerical  plane  wave  travel  over  many  computational  cells.  One 
modeling  technique  that  in  recent  years  has  become  one  of  the  most  popular  for  solving  Maxwell’s  equations  is  the 
Finite  Difference  Time  Domain  (FD-TD)  method  [l]-[3].  FD-TD  is  a  simple,  robust  and  easy  to  implement 
technique  for  solving  time-dependent  Maxwell  equations.  When  propagating  a  plane  wave  through  an  FD-TD  grid, 
there  are  three  commonly  used  algorithms:  the  total-field,  scattered-field,  and  separate-field  formulations.  Each 
technique  has  strengths  and  weaknesses  and  the  merits  of  each  are  discussed  in  [2-4].  All  techniques,  however,  have 
difficulty  with  plane  waves  traveling  over  many  cells.  In  this  paper,  a  new  scheme  in  conjunction  with  the  separate- 
field  formulation  is  derived  and  performance  results  are  given. 

II.  WAVE  NUMBER  ERROR 

All  numerical  solvers  have  inherent  errors  and  one  common  known  to  FD-TD  is  numerical  dispersion,  a  result  of 
phase  or  wave  number  error.  The  magnitude  of  the  wave  number  error  depends  upon  how  well  the  wavelength  is 
spatially  resolved  and  direction  of  travel  through  the  grid.  This  error  may  or  may  not  be  significant  depending  on 
the  electrical  size  of  the  mesh  and  the  number  of  wave  transits.  Recent  work  by  Nehrbass,  Jevtic  and  Lee  [5] 
discusses  a  technique  to  reduce  this  error  however,  another  related  error  arises  when  implementing  plane  waves 
which  is  the  focus  of  this  paper. 

A  plane  wave  by  definition  should  have  spatially  constant  field  values  within  all  planes  normal  to  the 
propagation  direction.  When  the  incident  wave  is  applied  to  the  grid,  a  3-D  numerical  plane  wave  is  excited  and 
propagates  through  the  grid.  If  at  each  frequency,  throughout  the  entire  bandwidth  the  applied  pulse  and  the 
numerical  wave  traveling  through  the  grid  do  not  have  identical  wave  numbers,  then  the  intended  plane  wave  will 
not  have  planes  of  constant  field  value.  This  difference  creates  fictitious  sources  along  the  total/scattered  field 
boundary.  In  many  cases,  the  error  grows  large  enough  to  corrupt  the  entire  grid  volume  yielding  inaccurate  results. 
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For  propagation  directions  parallel  to  one  of  the  grid  axes,  this  error  is  easy  to  rectify  by  using  interpolation  look-up 
data,  first  proposed  by  Tavflove  [3].  This  scheme  simply  interpolates  from  an  array  of  E-  and  H-field  values 
obtained  from  an  auxiliary  one-dimensional  FD-TD  grid.  Thus,  the  wave  applied  to  the  total/scattered  field 
boundary  can  be  exactly  replicated  when  propagated  parallel  to  one  of  the  axis.  In  this  paper,  this  technique  will  be 
referred  to  as  the  field  array  (FA)  scheme.  In  any  two  dimensional  diagonal  grid  direction,  the  dispersion  error  is 
small  making  it  simple  to  use  the  closed-form  incident  field  (CFIF)  scheme.  This  method  applies  the  field  values  to 
the  grid  via  an  analytical  equation.  However,  all  other  propagation  directions  are  not  as  easy  to  handle  with  either 
scheme  and  may  generate  large  errors  after  propagation  through  many  cells.  Some  recent  work  by  Oguz,  Gurel  and 
Arikan  [6]  greatly  reduces  the  error  by  filtering  out  high  frequency  noise  from  the  field  array  and  increasing  the 
spatial  resolution  of  the  FA  equation.  Their  work  shows  examples  of  three  dimensional  plane  waves  traveling  in  a 
two-dimensional  diagonal  direction  (kx=ky,  kz^^O)  and  primarily  focuses  on  a  single  frequency  excitation.  Instead, 
this  paper  develops  a  method  for  applying  a  wide  band  pulse  plane  wave.  The  FA  equation  is  reconstructed  with  a 
modified  cell  size,  permittivity  and  permeability  that  produces  a  close  approximation  to  the  3-D  numerical  wave 
number  band  over  the  entire  pulse  width.  For  the  rest  of  this  paper,  this  new  scheme  will  be  called  the  wide  band 
field  array  method  (WBFA).  Results  demonstrate  a  higher  fidelity  plane  wave  that  maintains  planes  of  near  constant 
field  values  normal  to  its  direction  of  propagation. 

III.  METHOD  DERIVATION 


This  section  first  briefly  discusses  the  dispersion  equation  and  then  derives  two  different  techniques  for 
optimizing  the  FA  equation  for  a  single  frequency  plane  wave  traveling  in  any  given  direction.  Each  method 
equalizes  the  wave  numbers  of  the  1-D  auxiliary  grid  and  the  3-D  grid  waves.  Then,  by  combining  these  two 
methods  via  an  iterative  algorithm,  the  FA  equation  is  further  optimized  to  simulate  a  band  of  wave  numbers  closely 
matching  those  of  the  3-D  numerically  propagating  wave.  For  a  given  direction,  the  wave  number  error  varies 
almost  linearly  with  the  number  of  cells  per  wavelength  as  long  as  the  shortest  wavelength  is  sufficiently  spatially 
resolved.  Thus  a  second  order  equation  can  be  used  to  match  a  band  of  frequencies. 

First  consider  the  transcendental  form  of  the  dispersion  equation  for  a  uniform  mesh  and  plane  wave  traveling  in 
a  given  direction  defined  by  0  and  (p: 


(Ar)'v^ 


2  COAt  2  ^  2  ^  2  ^ 


(3.1) 


where  Ax  =  Ay  =  Az  and 


For  a  given  frequency,  cell  size  and  time  step  satisfying  the  Courant  stability  limit,  there  is  a  resulting  numerical 
wave  number  that  satisfies  the  equation  above.  This  number  can  easily  be  found  by  an  iterative  guessing  method 
such  as  Newton’s  method.  Unfortunately,  the  numerical  wave  speed  a)/k  is  not  necessarily  equal  to  the  wave  speed 
as  it  ideally  should  be.  There  are  certain  conditions  for  which  they  are  almost  equal,  but  in  general  the  numerical 
wave  number  differs  from  the  exact,  analytical  k  that  will  produce  the  correct  wave  velocity  [2,  3]. 

Now  consider  the  1-D  dispersion  equation  that  applies  to  the  FA  equation  and  let  the  input  variables  be  the  cell 
size  and  wave  number  found  from  the  3-D  dispersion  equation  above: 
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Solving  for  the  speed  V  yields 
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(3.3) 


(Ax)^  sin^(-— ) 


(AO^sin^(^fc) 


e,P  and  ;Uj£,  of  the  FA  equation  can  be  adjusted  so  the  wave  number  matches  the  3-D  numerical  plane  wave 
number  found  in  equation  3.1.  Note  that  and  should  give  the  medium,  (usually  free  space)  wave 
impedance,  Z=  =4^xdIi^w  ■  By  implementing  the  new  and  ,  the  FA  equation  can  generate  the 

same  wave  number  as  the  3-D  numerical  plane  wave  for  a  single  frequency  yielding  a  much  improved  plane  wave. 

A  second  method  for  matching  the  wave  number  is  as  follows.  Let  the  input  variables  be  the  wave  number  and 
speed  found  in  equation  3.1  and  solve  for  the  cell  size  in  equation  3.2  using  an  iterative  method.  Again,  the  wave 
number  of  the  FA  equation  matches  that  of  the  3-D  numerical  plane  wave  for  a  single  frequency.  Both  methods  can 
apply  the  field  to  the  total/scattered  field  boundary  via  interpolation  from  the  FA  equation. 

Up  to  this  point,  the  FA  equation  wave  number  is  set  equal  to  the  3-D  wave  number  for  a  single  frequency  using 
one  of  the  two  methods  above.  However,  a  band  of  wave  numbers  needs  to  be  matched  to  accurately  propagate  a 
pulse.  This  task  is  a  little  more  complex.  Nonetheless,  by  combining  the  two  methods  above  it  is  possible  to  match 
an  entire  band  of  wave  numbers  with  a  high  degree  of  accuracy  by  employing  the  algorithm  shown  below.  The 
wave  number  at  the  upper  and  lower  frequencies  of  the  excitation  pulse  is  matched  by  adjusting  the  wave  speed  v,o 

and  the  cell  size  of  the  FA  equation.  The  scheme  is  as  follows: 


^min  =  Newton  iteration  of  3-D  dispersion  equation  ( Ax,  )  =  wave  number  at  lower  frequency  limit 

k  .  =  Newton  iteration  of  3-D  dispersion  equation  ( Ax,  )  =  wave  number  at  upper  frequency  limit 


Axjo  =  guess 


A  2  •  2  /  ^max  \ 

Sin  ( — - — ) 

2  „  (from  equation  3.3) 


k.  ,  n  =  Newton  iteration  of  1-D  dispersion  equation(  AXi^ ,  Vi^, ,  ) 


^min.lD  ^  ^min 


First,  the  wave  number  limits  are  found  from  the  3-D  dispersion  equation.  Since  FD-TD  commonly  uses  ten  cells 
per  wavelength  at  the  highest  frequency,  this  is  usually  a  good  upper  frequency  limit.  The  lower  frequency  can  be 
any  number  less  than  the  upper  limit  but  greater  than  zero  to  avoid  a  divide  by  zero  error  in  the  Newton  iteration. 
The  first  guess  for  the  cell  size  AXj^  is  the  3-D  grid  cell  size.  Next,  the  FA  equation  wave  speed  v,o  is  found 
directly  from  equation  3.3.  Now  that  the  upper  limit  wave  number  is  matched,  the  lower  limit  wave  number  is 
evaluated  via  Newton’s  iteration  method.  If  this  wave  number  is  not  equal  to  the  lower  limit  wave  number  of  the  3- 
D  dispersion  equation,  a  new  cell  size  is  approximated  and  the  loop  continues  until  both  the  upper  and  lower  wave 
numbers  match.  Successive  cell  sizes  can  be  approximated  by  interpolation  of  the  previous  two  guesses  and  the 
corresponding  wave  numbers  .  The  wave  number  variation  from  the  lower  to  upper  frequency  is  nearly 

linear  for  both  the  1-D  and  3-D  equations  thereby  producing  very  closely  matched  wave  numbers  and  phase  errors 
throughout  the  frequency  range.  Since  the  cell  size  and  time  step  of  the  WBFA  equation  differ  from  those  of  the  3- 
D  grid,  the  field  values  must  be  interpolated  precisely  correctly;  otherwise,  the  field  will  not  be  applied  accurately 
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but  instead  will  produce  large  errors.  A  nice  feature  of  this  method  is  the  efficiency.  Since  this  loop  is  executed 
during  preprocessing,  the  run  time  penalties  are  negligible. 

IV.  SINC  FUNCTION  EXCITATION  RESULTS 

To  measure  the  accuracy  of  the  WBFA  method,  a  series  of  plane  waves  were  launched  at  different  angles  of 
incidence.  In  order  to  attain  an  accurate  measurement  of  the  error,  an  empty  volume  with  no  scattering  objects  was 
used.  The  3-D  computational  domain  consisted  of  100  X  100  X  100  Yee  cells  encompassed  by  8-cell-thick 

perfectly  matched  layers  (PML)  with  a  normal  reflection  coefficient  R(0)  of  10"^  and  parabolic  conductivity 
profile.  The  cell  size  was  chosen  at  0.625cm  and  the  time  step  was  selected  to  be  just  under  the  Courant  stability 
limit  at  1 1.92  ps.  The  separate  field  formulation  was  employed  with  the  total  field  occupying  a  volume  of  88  X  88 
X  88  cells  and,  accordingly,  a  six  cell  thick  scattered  field  region  The  input  pulse  was  a  sine  function  with  a 
frequency  band  from  DC  to  4.8GHz.  Ten  cells  per  wavelength  resolved  the  highest  frequency.  The  peek  amplitude 
of  the  sine  wave  was  one  V/m  and  entered  the  volume  at  time  step  500.  All  tests  were  run  on  360MHz  SPARC 
Ultra  60’s. 

The  same  permittivity  and  permeability  values  for  the  3-D  FD-TD  equations  were  implemented  for  all  runs.  No 
optimization  for  one  propagation  direction  was  performed.  Instead,  an  optimized  permittivity  and  permeability  were 
computed  via  a  weighted  average  of  the  3-D  dispersion  equation  over  all  angles  of  incidence.  This  was  applied 
using  a  scheme  derived  by  Nehrbass,  Jevtic  and  Lee  which  greatly  decreases  the  average  wave  number  error  over  all 
directions  [5].  These  same  values  of  permittivity  and  permeability  were  used  in  the  3-D  dispersion  equation  (Eq. 
3.1). 

In  addition,  all  test  runs  were  performed  with  the  wave  number  of  the  FA  equation  matched  to  the  3-D  wave 
number  at  4.8GHz  by  adjusting  the  permittivity  and  permeability  of  the  FA  equation.  The  higher  the  frequency,  the 
lower  the  spatial  wave  number  resolution  and  the  greater  the  wave  number  error.  Thus,  the  wave  number  was 
matched  at  the  frequency  of  highest  wave  number  error. 

Ideally,  the  E-field  values  in  the  total  field  zone  should  be  exactly  equal  to  those  of  the  incident  pulse  of  the  FA 
equation  and  the  fields  in  the  scattered  field  zone  should  be  exactly  zero.  However,  there  will  be  error  due  to  the 
nature  of  finite  difference  numerical  methods.  The  results  shown  in  Figs.  1-4  give  the  maximum  E-field  error  over 
both  the  total  and  scattered  field  zones.  The  error  is  calculated  by  first  taking  the  difference  of  the  3-D  numerical  E- 
field  value  and  the  corresponding  E-field  value  on  the  I-D  auxiliary  equation.  Then,  this  difference  is  normalized  to 
the  maximum  amplitude  of  the  input  plane  wave  (I  V/m).  In  each  figure  the  test  results  for  the  FA  scheme  are  in 
part  (a)  and  the  test  results  for  the  WBFA  scheme  are  in  part  (b). 

The  first  test  run  launched  the  plane  wave  in  the  xy-direction  (6  =  90^  and  (p  =  45^ )  with  the  E  field  vector 
pointing  in  the  z-direction.  The  error  results  are  shown  in  Figure  1.  Initially,  the  error  rises  and  then  drops  after 
about  250  time  steps.  The  initial  excitation  of  the  wave  causes  high  frequency  noise  which  falls  off  after  it  passes. 
Then  the  error  rises  again  starting  at  about  500  time  steps  and  falls  at  about  750  time  steps.  At  500  time  steps  the 
peak  of  the  sine  wave  enters  the  grid.  This  peak  rises  far  above  the  previous  preceding  pulse  and  therefore  has  the 
same  effect  as  the  initial  excitation.  In  addition,  the  amplitude  is  larger  and  therefore  simply  has  a  larger  error.  As 
shown,  the  WBFA  scheme  results  in  an  error  about  three  orders  of  magnitude  less  than  the  FA  scheme. 

In  the  second  and  third  test  runs,  the  E  field  vector  has  nonzero  values  in  all  three  directions.  The  angle  of 

incidence  for  run  two  is  0  =  45^  and  (p  =  45^ .  In  this  case  the  E  field  vector  is  equal  in  magnitude  in  all  three 
components  and  the  peak  error  of  the  WBFA  method  is  about  two  orders  of  magnitude  below  the  FA  method  (Fig. 
2).  In  Figure  3,  the  angle  of  incidence  is  0  =  22.5^  and  (p  —  22.5^ .  In  this  case  the  E  field  vector  has  different 

magnitudes  in  all  three  components  and  the  peak  error  of  the  WBFA  method  is  more  than  an  order  of  magnitude 
below  that  of  the  FA  method. 

Figure  4  shows  the  results  of  the  final  run  with  an  angle  of  incidence  at  0  =  90^  and  ^  =  0^  .  As  explained 
before,  in  this  case  the  FA  equation  can  replicate  the  3-D  numerical  traveling  wave.  Here,  the  optimal  cell  size, 
permittivity  and  permeability  are  equal  respectively  to  those  of  the  3-D  grid  solver  and  the  FA  equation  yields  nearly 
identical  results.  The  resulting  error  is  due  to  the  interpolation  from  the  WBFA  or  FA  equation  as  it  is  applied  to  the 
3-D  grid. 

The  FA  and  WBFA  equation  accuracies  depend  on  the  propagation  direction.  Both  methods  have  highest 
accuracy  along  the  x-axis  and  decrease  in  accuracy  as  the  propagation  direction  is  changed.  However,  the  FA 
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equation’s  error  rapidly  grows  to  as  high  as  10  ’  (to  one  decimal  place  at  an  angle  of  0  -  90°  and  (p-45°.  On 
the  other  hand,  the  WBFA  method  attains  maximum  error  of  10  ^  at  an  angle  of  0  =  22.5  and  (p  =  225  . 


V,  CONCLUSIONS 

In  the  test  runs  performed,  the  WBFA  scheme  had  a  maximum  error  of  10“^  for  a  sine  pulse  containing  a 
frequency  band  from  zero  up  to  4.8GHz  at  ten  cells  per  wavelength.  For  many  applications  such  as  field 
distributions  in  building  environments,  this  is  more  than  acceptable  accuracy.  The  error  appears  to  be  greater  in 
directions  where  the  E  field  components  differ  in  magnitude  and  are  not  equal  to  zero.  Conversely,  the  FA  method 

has  errors  in  the  first  decimal  place.  „,r.T7A  *u  t- 

Accurate  propagation  of  wide  band  plane  waves  is  one  of  the  obvious  advantages  of  the  WBFA  method,  lime 
domain  methods  such  as  FD-TD  are  often  utilized  for  pulse  propagation.  Pulses  are  especially  convenient  for 
coupling  problems  since  the  input  pulse  and  the  coupling  points  of  interest  can  be  transformed  to  the  frequency 

domain  allowing  for  an  entire  frequency  band  of  coupling  results  with  one  code  run. 

Although  this  scheme  was  applied  to  the  separate-field  formulation,  the  same  WBFA  equation  could  be  applied 
to  the  scattered-field  formulation  as  well.  By  matching  the  wave  number  band  of  the  WBFA  equation  with  tha^t  of 
the  3-D  grid,  the  subtraction  noise  should  be  reduced  significantly.  This  is  especially  important  for  small  coupling 
problems. 
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